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ABSTRACT 



We present an improved method for determining the mass of neutron stars 
in eclipsing X-ray pulsar binaries and apply the method to six systems, namely 
Vela X-l, 4U 1538-52, SMC X-l, LMC X-4, Cen X-3, and Her X-l. In previous 
studies to determine neutron star mass, the X-ray eclipse duration has been 
approximated analytically by assuming the companion star is spherical with an 
effective Roche lobe radius. We use a numerical code based on Roche geometry 
with various optimizers to analyze the published data for these systems, which 
we supplement with new spectroscopic and photometric data for 4U 1538-52. 
This allows us to model the eclipse duration more accurately and thus calculate 
an improved value for the neutron star mass. The derived neutron star mass 
also depends on the assumed Roche lobe filling factor (3 of the companion star, 
where (3 = 1 indicates a completely filled Roche lobe. In previous work a range 
of (3 between 0.9 and 1.0 was usually adopted. We use optical ellipsoidal light- 
curve data to constrain (3. We find neutron star masses of 1.77 ± 0.08 M Q for 
Vela X-l, 0.87 ± 0.07 M & for 4U 1538-52 (eccentric orbit), 1.00 ± 0.10 M for 
4U 1538-52 (circular orbit), 1.04 ± 0.09 M Q for SMC X-l, 1.29 ± 0.05 M for 
LMC X-4, 1.49 ± 0.08 M for Cen X-3, and 1.07 ± 0.36 M for Her X-l. We 
discuss the limits of the approximations that were used to derive the earlier mass 
determinations, and we comment on the implications our new masses have for 
observationally refining the upper and lower bounds of the neutron star mass 
distribution. 

Subject headings: methods: numerical — pulsars: individual (Vela X-l, 4U 1538-42, 
SMC X-l, LMC X-4, Cen X-3, Her X-l) — stars: neutron — X-rays: binaries 
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1. Introduction 

A neutron star is a compact object that is the remnant of a massive star. The 
structure of a neutron star depends on the equation of state of nuclear matter under 
extreme conditions, specifically the relation between pressure and density in the neutron 
star interior. For a given equation of state, a mass-radius relation for the neutron star 
and a corresponding maximum mass can be derived. Many such theore tical equations o f 



state exist, ranging from "soft" — a mass upper limit as l ow as 1.5 (IBrown fc Bethe 



19941 ) — to "stiff" — a higher upper mass limit near 3 M Q (IKalogera fc Baym 



19961 ). The 



accurate measurement of neutron star masses is therefore important for our understanding 



of the equation o: 



by 



Kiziltan et al. 



state of matter in such high density situations (e.g., see the recent review 



20101 ). 



Eclipsing X-ray binary systems where the X-ray source is a pulsar can be ideal 
systems for a dynamical determination of the neutron star's mass. The orbital period, the 
semiamplitude of the optical star's radial velocity curve, the duration of the X-ray eclipse, 
and the projected semimajor axis of the pulsar's orbit (measured from the pulse arrival 
times) can be used to find the masses of both stars. We consider six systems where the 
required measurements have been made: Vela X-l, 4U 1538-52, SMC X-l, LMC X-4, Cen 
X-3, and Her X-l. The first five systems listed have OB supergiant companion stars, and 
the lattermost system has a somewhat less massive companion (~ 2 M ). 

In this paper, we present an improved method for determining the mass of neutron 
stars in eclipsing X-ray pulsars. In Section [21 we review the widely used analytic method 
for neutron star mass determination and reproduce previous results for the six systems. In 
Section [3j we present our numerical code based on Roche geometry to analyze the published 
data for the six systems. In Section HJ we present our results for each individual system 
with the incorporation of optical light curves. In Section [51 we discuss the implications of 
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our results for neutron star formation and the equation of state. 



2. Analytic Method 



2.1. Basic Equations 



and 



In this section we r eview the analytic method introduced by iRappaport Jossl (119831 ) 



Joss fc Rappaport 



(119841 ) that i s widely used to measure 



;he mass of the neutron star 



and its optical companion (e.g., see 



van Kerkwijk et al. 



1995a: 



van der Meer et al. 



2003). 



To begin, one can write the masses of the optical companion and the X-ray source (M opt 
and Mx, respectively) in terms of the mass functions: 

K X P(1 - e 2 f' 2 



M, 



opt 



2nG sin 3 % 



and 



M- 



x 



^o 3 pt^(l - < 



2\3/2 



i 

i + - 

Q 



(1) 



(2) 



2ttG sin 3 i 

where K x and K opt are the semiamplitudes of the respective radial velocity curves, P is 
the period of the orbit, % is the inclination of the orbital plane to the line of sight, e is the 
eccentricity of the orbit, and q is the mass ratio defined as 



q 



M n 



K, 



opt 



'opt Kx 

The values for Kx and P can be obtained very accurately from X-ray pulse timing 
measurements (the projected semi-major axis of the pulsar's orbit in light-seconds, ax sin i, 
is usually quoted in publications, from which one finds K x = 2nc a x sini/P), and optical 
and/or UV spectra can provide a value for K opt . 



(3) 



Assuming a spherical companion star, the inclination of the system is related to the 
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eclipse half-anglaj 9 e , the stellar radius R, and the orbital separation a by 



1 - (R/a) 



cos Of, 



(4) 



Following the approach in 



Rappaport fc Josa (119831 ). the radius of the companion star is 



some fraction of the effective Roche lobe radius 



R — /3Rl, 



(5) 



where Ri is the sphere-equivalent Roche lobe radius. We will refer to the fraction f3 as the 
"Roche lobe filling factor." Combining Equations fl4]) and (jSJ) yields 



smz 



1 - (3 2 (R L /a) 



cos 6L 



(6) 



Rappaport fc Josd (119831 ) provide an approximate expression for Rl/cl, the ratio of the 
effective Roche lobe radius and the orbital separation: 



R 

— A + B\ogq + C log 2 q, 
a 



(7) 



where the constants A, B, and C are 



A = 0.398 - 0.026ft 2 + 0.004ft 3 

B = -0.264 + 0.052ft 2 - 0.015ft 3 
C = -0.023 -0.005ft 2 . 



(9) 
(10) 



Here, ft is the ratio of the rotational frequency of the optical companion to the orbital 
frequency of the system. In other words, it is a measure of the degree of synchronous 
rotation, where ft = 1 is defined to be synchronous. These four expressions give the 



x The eclipse half-angle 9 e , or more specifically the semi-eclipse angle of the neutron star, 
represents half of the eclipse duration. 
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value of Rt, to an accura cy of about 2% over the ranges of < Q < 2 and 0.02 < q < 1 



( IJoss &: Rappaport 



1984 ). If the orbit is eccentric, the value of 8 is defined at 



periastron. 



19761 ). and at 



The star is assumed to have the same volume over its entire orbit (lAvnil 
any given phase outside periastron the value of 8 is adjusted accordingly. With a given 
eccentricity e, the orbital phase of the X-ray eclipse is determined by the argument of 
periastron to. 

For a given system, one can calculate the neutron star mass using Equations ([1]) through 
fTTU]) when given values of P, axsmi, 9 e , K opt , Q, and (3 (and if the orbit is eccentric, e and 
oj). It is possible to estimate fl by measuring the projecte d rotational velocity , v rnt s in i, 



van der Meer et al. 



fl2007h and 



of the optical companion star. This process is described in 
is employed there for the systems SMC X-l, LMC X-4, and Cen X-3. Finally, one must 
assume some value for the Roche lobe filling factor, /3. Many of the wind-fed systems 
discussed here are thought to b e close to filling their Ro c he lobes, and a range of 8 between 



0.9 and 1.0 is usually adopted ( IRappaport fc Joss 



1983 



van der Meer et al 



20071 ). Since 



all the measured input quantities are not known exactly, a simple Monte Carlo technique 
may be used to derive the most likely values of M x , M opt , and i, and the corresponding la 
confidence limits. 

To determine neutron star masses in this manner for the six eclipsing systems where 
all of the necessary quantities are known or estimated (see Tables [1] and [2] and references 
therein), we assume that any value within the range 0.9 < 8 < 1.0 is equally likely. We then 
use a Monte Carlo technique to generate a distribution of resulting neutron star masses, 
as shown in Figure [TJ This provides a clear visual representation of uncertainties through 
the shapes of each distribution, and it provides an estimate of the "most likely" mass for 
each system. The neutron star masses and system inclinations found in this manner are 
presented in Table [31 The masses derived for SMC X-l, LMC X-4, and Cen X-3 agree very 
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well with those cited in Ivan der Meer et al.l (120071 ) 



We note that there is no physical solution for either Vela X-l or 4U 1538-52 when using 
this technique with an eccentric orbit, because the quantity in Equation ([6]) is larger than 
unity. Adjusting the argument of periastron u and/or the eclipse duration 9 e can force a 
solution. However, for 4U 1538-52, no solution exists within the la uncertainties of u and 
9 e . This discrepancy arises due a high inclination and is discussed further in Sections 14.11 



and 14.21 As a workaround, when employing this ana. 



width 



33 ± 3°, from Ivan Kerkwijk et al. 



orbit (e = 0) as in 



ytic technique we use a larger eclipse 



1995bl ) for Vela X-l and adopt a circular 



Clark! ( 12000( 1 for 4U 1538-52. 



2.2. An Examination of the Approximations 



The analytic method presented in Section 12.11 is straightforward, and is easy to 
implement on a computer. However, this method relies on two approximations: 



i. The computation of the effective Roche lobe radius, Rl/cl, from Equations (J7])- (|10p . 
and 

ii. The computation of the X-ray eclipse duration, 26 e , from Equation (J6]). 



We use the Eclipsing Light Curve (ELC) code of lOrosz fc Hauschildtl ( 12 00 01 ). which is based 
on Roche geometry, to test these two approximations, and we discuss each one in turn. 

Strictly speaking, "Roche geometry" applies only to binary systems with circular 
orbits and co-rotating stars. Numerous authors have presented generalizations of the 
Roche pot ential to account for s i tuations in which one or both of these assumpt ion are not 



met (e.g., 



Avni fc Bahcall 



1975 



Avni 



1976 



Wilson 



1979 



Avni fe Schiller 



1982). These 



generalizations do not fully describe complete dynamics of the star, and as a result some 
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small approximations are involved (e.g. Wilson 1979). The main assumption is that the 
timescale for the internal motions of the star that are required for the star to adjust to 
the varying potential is considerably shorter than the orbital period. Given this, one can 
compute an ef fective potential locally at each orbital phase without significant inconsistency 



((Wilson 



19791 ) . Although these generalized potentials are widely used, it is not known how 
well they work in practice. For most of the systems discussed here, the orbits are circular 
and the stars rotate close to the sync hronous rate, so the modified potential we use (the 



ELC code uses the potentia 
note that 



given m 



Wilson 



19791 ) should be fairly accurate. Finally, we 



Rappaport & Jossl (119831 ) also adopted a modified "Roche potential" in their 
analysis — their fitting functions are approximations to numerical integrations of the critical 
potential surface. Therefore, any systematic error ELC would have owi ng to improper 



genera lizations of the Roche model would also be present i n the work of 



Rappaport fc Joss 



( 119831 ) and others that used these fitting functions such as Ivan der Meer et al. 



(120071 ) 



The shape and size of the critical potential surface (hereafter the "Roche lobe" ) depend 
only on the mass ratio q and the parameter Q. When given q and it is straightforward to 
define the equipotential surface (from the value of the gravitational potential at the inner 
Lagrangian point) and to numerically integrate its volume. The sphere-equivalent radius 
Rl then follows. We define a large grid of points in the q-fl plane and compute values of 
Rl, and compare them with t he values of Rr, fou n d from Equations (jTJ) through (flOl) . The 



results are shown in Figure [2j 



Rappaport fc Jossl (119831 ) claim an accuracy of their fitting 



functions of about 2% over the stated range (0 < Q < 2 and 0.02 < q < 1), and our results 
confirm this. 

From Equation (Q, one can see that the duration of the X-ray eclipse depends on the 
inclination i, the Roche lobe filling factor /3, and Ri/a, which is a function of the mass ratio 
q and the parameter Q. The ELC code can be used to compute the duration of the X-ray 
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eclipse for a given geometry. Rather than u sing ray tracing to determine whether a point 



is eclipsed by the companion star (e.g., see 



Chanan. Middleditch. fc Nelson 



19761 ). ELC 



locates the limb of the star to high accuracy by testing the viewing angles of each surface 
element. ELC then uses bisection to find, at each latitude row, the longitude of the point 
that has a viewing angle of fi = 0. Once found, the points on the limb define a polygon in 
sky coordinates. At that same phase, the location of the X-ray source (assumed to be a 
point source) in sky coordinates is determined. A simple test is used to determine if the sky 
coordinate of the X-ray source is inside or outside the polygon defined by the horizon of 
the star. ELC uses another bisection routine to find the orbital phase of the X-ray eclipse 
ingress to high accuracy. If the orbit is circular, the eclipse half angle 8 e is equal to the 
ingress phase. If the orbit is eccentric, the X-ray eclipse egress phase is also computed, and 
the eclipse half angle is computed from both the ingress and egress phases. 

Setting ft — 1, we use ELC to compute the full duration of the X-ray eclipse for a wide 
range of values in the q-i plane, using (3 = 1.0 and [3 = 0.9 and assuming a circular orbit. 
The full eclipse duration was also computed from Equation (jHJ), and the difference between 
the numerically computed value of 28 e and the analytically computed value of 28 e was 
determined. Figure [3] shows the differences for (3 = 1 and Figure H] shows the differences for 
(3 = 0.9. The differences can be quite extreme: They are in excess of 10° for small mass 
ratios and large inclinations and are less than —10° for grazing eclipses. 

The approximate locations in the q-i plane for three systems (SMC X-l, LMC X-4, and 
Cen X-3) are also shown in Figures [3] and HI Since Her X-l is a low mass X-ray binary, its 
mass ratio does not appear within the limits of the two figures. Vela X-l and 4U 1538-52 
are excluded because of their eccentric orbits. The instantaneous Roche lobe filling factor 
of an eccentric system during X-ray eclipse will be less than the value of (3 as calculated 
for non-eccentric systems. When — 1, all three systems shown in Figures [3] and H] are 
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near the contour denoting zero difference. When f3 = 0.9, all three systems are above the 
zero-difference contour, which indicates that the eclipse durations computed numerically 
are longer than those computed analytically. 

To illustrate the differences between the analytic and ELC results, Figures [5] and [6] 
show a system resembling Cen X-3 in sky coordinates, where the Roche lobe filling factor 
is (5 = 0.9. Figure |5] demonstrates that the companion star is not spherical. The parts of 
the star near its equator extend beyond the circle denoting the volume-equivalent sphere. 
Hence one would expect that the duration of an X-ray eclipse for very high inclinations 
would be longer than what one would compute from the analytic approximations, and a 
glance at Figures [3] and H] confirms this. In a similar manner, Figures [3] and H] show that for 
grazing eclipses (i.e., lower inclinations), the numerically computed durations are shorter 
than the analytic approximations. One can see from Figure [5] that the polar regions of the 
companion star are inside the circle denoting the volume-equivalent sphere, and as a result 
there would be inclinations at which the analytic approximations indicate X-ray eclipses 
when in fact none occur. 

Figure [6] shows a magnified view of Figure [5] near the neutron star eclipse. In this 
example, the egress phase of the X-ray eclipse is very close to 33°, since the neutron star 
is just crossing the limb of the companion star. The limb of the analytically computed 
volume-equivalent sphere is well inside the limb of the star. The egress phase of X-ray 
eclipse would be near 32° when the analytic expressions are used, since the neutron star is 
just crossing the limb of the sphere at that phase. Thus, in this example, the full duration 
of the X-ray eclipse computed numerically is a full 2° longer than the duration computed 
analytically. 

Having tested both approximations, we conclude that the ELC code does offer a 
significantly more accurate representation of the physical system than the analytic method 
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presented in Section 12.11 The size of the effective Roche lobe radius as a function of q and 
Q is relatively well represented by the analytic formulae and ELC offers only a modest 
improvement. However, the difference in X-ray eclipse duration can be quite extreme (as 
much as ±10°) and has a direct effect on the calculation of the neutron star mass. 

3. Numerical Method 

We use the ELC code and its various optimizers to analyze the data given in Table 
1 and optical light curves (see Section H]) to derive the neutron star masses and their 
uncertainties. The ELC code has two advantages here. First, Roche geometry is used (i.e., 
no approximations are used to find the effective Roche lobe radius or the X-ray eclipse 
duration, as discussed in Section [272]) . Second, when using ELC, one can make use of any 
number of other sources of information about the system, such as optical light curves. 
When the geometry is specified, one can compute various observable properties of the 
system and compare them with the observed values using a \ 2 or similar test. One can 
then find the family of geometries that best match the observed quantities, and from those 
geometries the masses and other system parameters follow. 

To begin, the orbital period P (which is known to high accuracy) and the orbital 
separation a give the total mass of the binary via Kepler's Third Law. Specifying the mass 
ratio q then gives the component masses. The shape of the companion star is determined 
when the Roche lobe filling factor and the parameter Q are given. The shape of the 
orbit, if eccentric, is determined from the eccentricity e and the argument of periastron 
oj. Finally, when the inclination i is given, it is possible to find the K- velocities of the 
components, the rotational velocity of the companion star, and the duration of the X-ray 
eclipse. Thus we initially have an eight-dimensional parameter space (P, a, q, (3, Q, e, u, i) 
to search. However, the search of the parameter space can be simplified. First, we assume 
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the period P is known exactly and fix it at the appropriate value for each system. Likewise, 
the value of the projected semimajor axis of the pulsar's orbit ax sini is usually known to 
high accuracy by measuring the X-ray pulse arrival times, and is also held fixed. Next, the 
parameters a and q can be computed from ax sin i and the K- velocity of the companion 
star: 

2ttc(1 -e 2 ) 3/2 a x sini 
K x = p (11) 

q = §Hi (12) 
&x 

a = ^1 H — j c ax- (13) 

These simplifications give us a six-dimensional parameter space (K opt , (3, fl, e, u, i) to 
search. 

ELC's various optimizers allow the user to specify a wide range of values for each 
parameter as well as sets of data points related to the physical system, such as light curves 
or radial velocity curves. ELC forms random sets of parameters and uses each set to 
compute a model. The "fitness" of each model is defined using a \ 2 merit function: 



2 f fl e (mod)-fl e (obs) y 

O0e / 



X 

+ 



+ 



+ 



v rot sin i (mod) — v rot sim(obs) 



K opt (mod) - K opt (obs) \ 2 



2 



a K opt 



1 



e mo 



d) — e(obs) V 



(mod) — w(obs) \ ' 



+ , ^ , _ (14) 



Here, the notation (mod) means the quantity computed from the model, the notation (obs) 
means the observed quantities, and the notation ctq indicates the la uncertainty of the 
observed quantity. We use the values of v rot sini and oj given in Table [2j For Her X-l, we 



13 



do not constrain v ro t sin i but assume synchronous rotation. It might seem strange to have 
the terms involving K opt , e, and u in the merit function since they are input parameters for 
the model, but our experience has been that the optimizers perform better when the input 
parameters are drawn from a uniform distribution. 

Once the fitness of a given model is determined, new parameter sets are constructed 
using either a Monte Carlo Markov chain optimizer or a genetic algorithm optimizer. In 
the latter case, a "breeding" technique is used that is based on the "survival of the fittest" 



( lOrosz et al. 



20021 ). The probability of "breeding" is based on a model's fitness. Random 
variations (i.e., "mutations") are introduced into a small fraction of the breeding events, and 
the process of br eeding a new populat ion and evaluating its members is repeated over many 



generations. See ICharbonneaul (jl995f ) for a more detailed discussion of genetic algorithms. 
In both the Monte Carlo Markov chain and the genetic algorithm, the fitness of each new 
parameter set is determined and the process is repeated until convergence is achieved. 

Ultimately, we compute hundreds of thousands of models. Each model has an 
associated value of x 2 an d various derived parameters including the component masses, 
system inclination, etc. A lower x 2 value indicates higher fitness. In this particular case, the 
minimum possible value of x 2 is zero since the there are more linearly independent input 
parameters than observed quantities. To define the la limits, the family of models where 
X 2 < 1 is found and the distributions of the various parameters of interest are constructed. 

We note that the value of ax sin % given for 4U 1538-52 has a relatively large uncertainty. 
To account for this, we modified the genetic code to allow for a range of ax sin i values 
drawn from the appropriate Gaussian distribution to be used. Although this modification 
hardly made a difference in the output parameter distributions, our results do account for 
the uncertainty in ax sini. 

We performed a preliminary analysis of the six systems for a range of 0.75 < /3 < 1 and 
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those results are shown in Figured The system inclination is inversely correlated with ft, 
and the lower bound on ft for each system corresponds roughly to i — 90°. Hence not all 
systems can have masses for all of the values of ft considered. Two things are apparent in 
Figure [7J First, the neutron star mass is highly dependent on the choice of ft, and second, 
the numerical and analytic results can differ in opposite senses to varying degrees depending 
on the choice of ft. 



To improve upon the technique described in Section [3j we consider optical light curves 
for five systems in addition to the parameter values in Tables [H and EJ The shape of the 
ellipsoidal light variations from the companion star depends on the inclination i, the mass 
ratio q, the parameter Q, and the companion star's Roche lobe filling factor ft. Since the 
first three parameters are already well-determined from the width of the X-ray eclipse and 
the K- velocities of each component, ft should be quite well constrained by including optical 
light curves in the analysis. Such an analysis, which we now present, is trivial to do using 
ELC and extremely difficult to do using the analytic approximations. 

Adding new observables to the ELC analysis adds terms to the x 2 merit function as 
originally shown in Equation (I14p. We effectively have 



where the final term incorporates a set of iV observations defined by observable quantities yi 
(e.g., an optical light curve with N data points). Similar terms may be added for additional 
sets of observations (e.g., a radial velocity curve). 

All of the folded optical light curves discussed below are presented in Figures [8] and [9] 
with the best fit model for each system from ELC. For the five systems we analyze in this 



4. Optical Light Curves 




(15) 
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manner, the best fit solution includes an accretion disk around the neutron star. Adding a 
disk increases the depth of the secondary eclipse and results in a better fit in all cases. For 
comparison, we have plotted both the best fit model and the same model with the light 
from the accretion disk subtracted in Figures [8] and [H Our final neutron star masses are 
presented in Table H] and Figure [10] with the corresponding analytic masses. Rather than 
using a range of 0.9 < /3 < 1 for the analytic cases, we use the value for /3 returned by the 
best fit ELC model for consistency. 



4.1. Vela X-l 



Vela X-l has a n 8.96 day orbital period and an eccentric orbit with e = 0.0898 ± 0.0012 
LJ) . The argument of periastron for the companion star is u = 332.59±0.92° 



Barziv et al. 


2001 




1 

(Bildsten et al.lll997n 



(IBildsten et al.l 119971 ). where u) = ujx + 180. The optical V light curve shown in the upper 



l eft panel o 



(jPojmanskv 



Figure E] for Vela X-l is binned data from the All Sky Automated Survey 



20021 ) . In our preliminary analysis of Vela X-l (e.g., Figure [TP , we adopt a 



semiduration of the X-ray e clipse 6 e = 33 ± 3° ( 



the more precise value from 



Kreykenbohm et al. 



van Kerkwijk et al. 



1995bJ). We later use 



(120081 ) . 6 e = 34.135 ± 0.5°. Unfortunately, 



the analytic technique from Section [2] does not arrive at a physical solution with this longer 
eclipse duration due to the system's high inclination. Specifically, solving for the inclination 
i is not possible using Equation ([6]) as sini > 1 for all the Monte Carlo simulations. 
The numerical ELC code does not have this same limitation. As a workaround, we keep 
9 e = 33 ± 3° for all instances of the analytic case instead. Our final derived mass for Vela 
X-l is 1.77 ± 0.08 M with a system inclination of 77.8 ± 1.2°. 
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4.2. 4U 1538-52 



4U 15 38-52 has a 3.73 day orb ital period and most likely an eccentric orbit. I Clark 



( 12000( 1 and 



Mukherjee et al 



whi ch are in agreement. Howeve r, 



(12007/) give e = 0.174 ± 0.0 



Makishima et al. 



and 



van Kerkwijk et al. 



5 and 0.18 ± 0.01, respectively, 



(119871 ) give a mu ch low e r e = 0.08 ±0.05 



Clark! (120001 ) provides 



(Il995bl ) adopt e = in their analysis; even 
an alternate set of fit parameters for a circular orbit. As with Vela X-l, we find no physical 
analytic solution for 4U 1538-52 with an eccentric orbit and its reported eclipse duration 
(9 P . = 28.5 ± 1.5°) due to the system's high inclination. Therefore, we follow the approach of 



Clark! (120001 ) and present numerical results for both an eccentric orbit and a circular orbit. 



In addition, 



Clark! (I2000h and 



Mukherjee et al. 



( 120071 ) give significantly different 



values for u, the argument of periastro n: 244 ± 9° and 22 ± 12° for the optical source 



respectively. Following the approach of 



Zhang et al. 



(120051 ). we have extrapol ated a 



decr ease of oj based on the time difference between the observations made by 



and 



Mukherjee et al. 



Clark! (12000h 



mear 



(120071 ). We estimate u is decreasing by 0.0010 ± 0.0006° per day and 
that our subsequent observations should therefore have u = 198 ± 14°. This is the value we 
adopt in the numerical model. Even though using a lower u nudges the parameter space 
closer to a single physical analytic solution, the reported eclipse duration is still too short 
to allow sini < 1 (as in Equation [6]) , and we are forced to retain a model where e = for 
the analytic case. 

Due to the above disc repancies, the reportedly low neutron star mass (~ 1M C 



(Ivan Kerkwijk et al. 



1995b! ) . and lack of a published light curve, we found this system 



especially worthy of our attention. 

Optical light curves in BVI for 4U 1538-52 were obtained at the Cerro Tololo 
Inter- American Observatory on the 1.3 m SMARTS telescope with the ANDICAM in 
June - September 2009. There are a total of 39 images in each filter taken on different 
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nights, each with a 60-second exposure time. Standard pipeline reductions were done on 
all images, and differential photometry was performed in IRAF for the target and several 
comparison stars. Light curves are shown in Figure O and observations in all three filters 
were incorporated into our numerical analysis. 

Twenty-one high resolution spectra of 4U 1538-52 were also taken on several nights in 
July and August 2009 at Las Campanas Observatory on the 6.5 m Clay Magellan telescope 
with the MIKE spectrograph. Standard pipeline reductions were done on all frames, 
including heliocentric corrections. The spectroscopic analysis was performed in IRAF using 
cross-correlation of the blue half of the spectrum (4750 - 4950 A) with a model BO star. 
Radial velocities were then computed and incorporated into our numerical analysis. The 
radial velocity curve is plotted in the bottom panels of Figure |9] with the adopted ELC 
model solution. 

We ran two sets of model fits. We assumed an eccentric orbit with eccentricities in a 
narrow range around e = 0.18 and also a circular orbit. We find K opt = 14.1 ±1.1 km s _1 for 
an eccentric orbit and K opt = 21.8 ± 3.8 km s _1 for a circular orbit. The difference between 
these two measurements is due in part to the noi s y dat a and incomplete phase coverage 



near phase 0.25. For comparison, iReynolds et al. 



{19921) found K opt = 19.2 ± 1.2 km s" 1 



(uncorrected) and K opt = 19.8 ± 1.1 km s _1 (corrected for tidal distortions) assuming a 
circular orbit. Our final mass for 4U 1538-52 is quite low: 0.87 ± 0.07 M & for an eccentric 
orbit and 1.00 ± 0.10 M & for a circular orbit. If the orbit is indeed eccentric, this is an 
extremely low neutron star mass. 
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4.3. SMC X-l 



SMC X-l has a 3.89 day orbital period. It also exhibits a superorbital X-ray cycle 



that is probably caused 



Priedhorsky &: Holt 



1987) 



Dy precession of either the accretion disk or neutron star (e.g., 
The optical V light curve fo r SMC X-l shown in the lower left 



panel of Figure [8] is from 



van Paradijs fc Kuiperl (jl984|). Our final derived mass for SMC 



X-l is 1.04 ± 0.09 M & with a system inclination of 68.5 ± 5.2°. We discuss the implications 
of this low mass result and that of 4U 1538-52 in Section 



4.4. LMC X-4 



LMC X- 4 has a 1.41 day orbital per i od and , like SMC X-l, also exhibits a superorbital 



X-ray cycle. 



Heemskerk fc van Paradijd (119891 ) performed an extensive study of the 



long-term variations in X-ray flux of LMC X-4 and concluded it contained a warped, 
precessing accretion disk. The optical B light curve for LMC X-4 shown in the upper 



right panel of Figure [8] is from 



Ilovaisky et al 



( 119841 ) . Since no data table was available, 



the points were extracted from their "X-ray OFF states" plot using DEXTER available 
via the SAO/NASA Astrophysics Data System. Our final derived mass for LMC X-4 is 
1.29 ± 0.05 M Q with a system inclination of 67.0 ± 1.9°. 



4.5. Cen X-3 



Cen X-3 has a 2.09 day orbital perio d. The optical V 



the lower right panel of Figure |8] is from 



van Paradijs et al 



i ght cu rve for Cen X-3 shown in 



(119831 ). Our final derived mass 



for Cen X-3 is 1.49 ± 0.08 M with a system inclination of 66.7 ± 2.4°. 
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4.6. Her X-l 



Her X-l has a 1.7 day orbital period and a much lower companion star mass than 
the other five systems of interest (~ 2M odot ). Like SMC X-l and LMC X-4, Her X-l also 
exhibits superorbital X-ray cycles. 

No optical light curve was used for Her X-l because the optical iT-velocity has a very 
large uncertainty (90 ± 20 km s _1 ) due to uneven X-ray heating of the companion star. As 
a result, any spectral lines measured do not on average emanate from the star's center of 
mass. An optical light curve would therefore be ineffective in constraining the fit. Further, 
the companion's Roche lobe filling factor is likely very clos e to 1 (between 0.95 and 1) as 



adopted by and discussed in 



van Kerkwijk et al. 



fll995bf ) and 



Bahcall fc Chester! ( 119770 . We 



therefore set (3 = 1 in our analysis. Our final derived mass for Her X-l is 1.07 ± 0.36 M & 
with a system inclination greater than 85.9°. 



5. Discussion 

We present our final derived neutron masses for all six systems in Figure [TO] These 
values are also given in Table HI In every case where we have incorporated light curves into 
the analysis, our la error bars are smaller than those of the analytically-derived masses. It 
is clear from this figure that these six neutron stars have different masses that span a range 
from as low as 0.9 M Q to as high as 1.8 M & . 



Kiziltan et al 



(l2010f ) recently provided a comprehensive review of mass determination 

for neutron stars in double neutron star binaries and in neutron star-white dwarf (NS-WD) 

binaries. They did not include neutron stars in X-ray binaries since the mass determinations 

generally have larger uncertainties than the typical uncertai nties for the double ne utron star 

I ~- I J I 

and NS-WD systems. The lowest mass neutron star in the iKiziltan et al.l (120101 ) sample 
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with a small uncerta inty is the companion to PSR J1756-2251 with M = 1.18 ± 0.03 M 



( IFaulkner et al. 



20051 ). The largest mass wi th a small uncertainty is PSR J1614-2230 



M = 1.97 ± 0.04 M Q ( IDemorest et al. 



2010J). Other high mass ne utron stars include 



PSR B1516+02B with M = 2.10 ± 0.19 ( Kiziltan et al. 



with M = 1.91±q°1,Mq (iKiziltan et al 



20101 ). PSR J1748-2446I 



2010). and possibly the "Blac k Widow Pulsar' 



PSR B1957+20 with M = 2 AO ± 0.12, M (Ivan Kerkwiik et al. 



20101). In the case of 



the lattermost system, the systematic errors are potentially large owing to the extreme 
irradiation suffered by the pulsar's evaporating companion. The large and secure (low- 
mass of PSR J1614-2230 rules out many of the so-called soft equations of state 



uncertainty 


mass 


(Ozel et al. 


2010) 



As discussed by 



Kiziltan et al. 



(120101 ). the core mass of a star needs to exceed 



the Chandrasekhar mass if it is to end up as a neutron st ar. The exact value of the 



Chandrasekhar mass depends on the electron fraction, and 



Kiziltan et al 



( I2010h gives a 



plausible range of possible neutron s tar birth masses of 1.08 < Mbi rt h ^ 1.57 M Q . Finally, 



as discussed by Kiziltan et al. 



(120101 ) . the millisecond pulsars in the NS-WD binaries must 
have accreted some mass in order to end up with millisecond spin periods. They give 
a range of 0.10 < M acc < 0.20 M based on angular momentum considerations and on 
plausible mass transfer rates in the X-ray binary phase. They conclude that neutron stars 
with masses less than about 1.1 M & would be unusual since it would be difficult to exceed 
the Chandrasekhar mass, while neutron stars with masses above about 1.8 Mq must have 
had either a prolonged stage of mass transfer at an unusually high rate or an unusually high 
mass when they formed. 

Until recently, the neutron star in Vela X-l has been at the high end of neutron star 
mass measurements. Unfortunately, the systematic errors in the mass determination are 
difficult to minimize owing to the non-radial pulsations in the B-giant companion. Using 
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the K- velocity for Vela X-l given in iBarziv et al.l ( 1200 ll ). the mass of its neutron star is 



1.77 



0.08 M and the companion star fills its Roche lobe at periastron. 



Quaintrell et al 



( 120031 ) derive a slightly 



21.7 ± 1.6 km s 1 from 



higher K-ve 



Barziv et al. 



ocity for Vela X-l (22.6 ± 1.5 km s 1 compared to 



20011 ). and consequently they derive a higher mass for 
the neutron star using the analytic approximations (2.27 ± 0.17 M for (3 = 1). When we 
used this higher i^-velocity in the ELC code, the best fit gave a more conservative neutron 
star mass of 1.84 ± 0.06 M . However, the overall \ 2 was a bit worse (x 2 = 30.96 compared 
to 29.52), and we note that this mass does barely fall within the la uncertainty of our 
adopted value. In spite of the uncertainties, it appears likely that the neutron star in Vela 
X-l has a relatively high mass that is comparable to the masses found for PSR B1516+02B, 
PSR J1748-2446I, and the Black Widow Pulsar. It is interesting to note that the companion 
star in Vela X-l is high mass and will most likely also produce a neutron star. If so, then 
Vela X-l cannot end up in a state similar to the Black Widow Pulsar or to the NS-WD 
binaries PSR B1516+02B and PSR J1748-2446I. 

On the lower mass end, we have 4U 1538-52 at 0.87 ± 0.07 M (eccentric) or 
1.00 ± 0.10 M (circular) and SMC X-l at 1.04 ± 0.09 M . Further spectroscopic 
observations are needed to better establish the shape of the orbit in 4U 1538-52 and to 
better establish the K- velocity. In the case of SMC X-l, the observational uncertainties are 
much smaller. Its mass is about 1.5a smaller than the companion to PSR J1756-2251, which, 
as noted above, is the least massive neutron star known with a secure mass determination. 
The pulse timing properties of SMC X-l are very well known, and the only other main 
quantities that could potentially change are the i^-velocity of the companion star and the 
duration of the X-ray eclipse. To see what has to change in order to have a neutron star 
mass of ~ 1.2 M & in SMC X-l, we performed a series of fits to the light curve where the 
i^-velocity and eclipse duration were fixed at various values. To drive the neutron star 
mass to higher than 1.2 M , the optical .ff -velocity needs to be higher (K ovt <i 22 km s _1 or 
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almost 2cr higher than observed) and the eclipse duration needs to be shorter (8 e < 27.5° or 
about la smaller than observed). These simulations are presented in Table [5J If confirmed, 
the low masses for the neutron stars in SMC X-l and 4U 1538-52 would indeed challenge 
neutron star formation models. 

There are a few ways to modestly improve the accuracy of our neutron star mass 
determinations. Since the mass of these neutron stars is proportional to K^ pt , improvements 
in the measured velocity curves will improve the mass determinations. For three of the 
systems (SMC X-l, LMC X-4, and Cen X-3), this would involve the acquisition and analysis 
of a significant number of additional high resolution optical spectra, work which we did 
perform in the case of 4U 1538-52 (see Section I4~2l) . Improvements in the K- velocity for 



Vela X-l wi 



companion ( jBarziv et al 



1 be difficult owing to the presence of no n-radial pulsations in its B-giant 



2001 



Quaintrell et al 



20031 ) . Likewise, improvements in the 



K- velocity of Her X-l will be difficult owing to the relatively strong X-ray he ating of its 



low- mass companion (e.g., 



Crampton &: Hutchings 



1974 



Reynolds et al. 



19971 ). Modest 



improvements in the measured eclipse widths and X-ray timing properties for many of the 
systems could be made by observing the sources with greater time coverage. 



MLR and JAO gratefully acknowledge the support of NSF grant AST-0808145. We 
also thank the anonymous referee who provided constructive feedback and encouraged us 
to include the full optical light curve analysis in this work. 
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Fig. 1. — Resulting probability distributions (histograms) from Monte Carlo simulations 
using the analytic method for all six systems as a function of neutron star mass. We assume 
that any filling factor 0.9 < /3 < 1.0 is equally likely. The mean value ±lcr is given above 
each peak. All input parameters for these distributions are given in Tables U and [21 We not e 



that our mass for Cen X-3 agree very well with that found by Ivan der Meer et al. 



(200 



an d the masses for SMC X- l and LMC X-4 also agree exactly when the values of 9 e given 



in 



van der Meer et al. 



(120071 ) are used. We assume e = for 4U 1538-52. 
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Effective radius difference (%), numerical-analytic, /S=1.0 




Mass Ratio (M x /M opt ) 

Fig. 2. — Percent difference between the numerically computed effective Roche lobe radius 
and the analytically computed one (Equations (T7|)- (fTUl) ) as a function of the mass ratio q 
and the parameter Q. The largest deviations are about 1.5% for mass ratios near 0.02 and 
1.0. 
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Eclipse duration difference (deg), numerical— analytic, /S=1.0, 0=1.0 
a> I ' ' ' ' 1 ' 1 ' ' 1 ' 1 ' 1 1 1 ' ' 1 1 ' 1 1 — 
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03 




0.05 0.1 0.15 0.2 0.25 



Mass Ratio (M x /M t ) 

Fig. 3. — Difference (in degrees) of the full X-ray eclipse duration computed numerically and 
the full X-ray eclipse duration computed analytically via Equations (16T)- f[HJj) as a function 
of the mass ratio q and the system inclination i. Here we set the parameter Q = 1 and the 
Roche lobe filling factor f3 = 1. The black curve denotes the contour corresponding to a 
difference of zero degrees. The locations for three systems are shown: C=Cen X-3, S=SMC 
X-l, and L=LMC X-l. 
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Cen X-3, j8 = 0.9, i=76.35°, orbital phase 34 




-10 -5 5 10 



x (R Q ) 

Fig. 5. — Representation in sky coordinates of a system resembling Cen X-3. The Roche 
lobe filling factor is /3 = 0.9, the inclination is i — 76.35°, and the orbital phase is <fi = 34°. 
The rotation axis of the companion star and the angular momentum vector of the orbit 
are parallel to the y-axis. The red circle denotes a sphere with the same volume as the 
companion star. The dot near the "10 o'clock" position marks where the neutron star is 
located at this phase. 
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Cen X-3, (8=0.9, i=76.35 




-10.5 -10 -9.5 



x (R ) 

Fig. 6. — Magnified view of the schematic diagram shown in Figure |5j The solid red line 
denotes a sphere with the same volume as the star, computed numerically. The dashed 
line denotes a volume-equivalent sphere computed using the analytic method via Equations 
fl6|)- f|T0|) . The location of the neutron star is shown at four orbital phases. In this example, 
the true egress phase from X-ray eclipse would be very close to 33° while the egress phase 
computed using the analytic approximations would be very close to 32°. 
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Fig. 7. — Preliminary neutron star mass versus the companion Roche lobe filling factor (3 
for all six systems. Analytic indicates masses derived using approximations (see Section [2]) 
while Numerical indicates masses derived using the ELC code based on Roche geometry (see 
Section Optical light curves have not been included in this analysis. Each numerical 
solution is equally likely in that the eclipse width in the model matches the observed eclipse 
width exactly. All residuals are in the sense Numerical-Analytic. For simplicity, 4U 1538-52 
is approximated as a circular system (e = 0). 
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Orbital Phase 



Fig. 8. — Optical light curves for four systems. The solid lines represent the best fit ELC 
model for each system, all of which include an accretion disk around the neutron star. 
The dotted lines depict the models with the light from the accretion disk subtracted for 



comparison. The 



(jPojmansky 



2002 



ight curves for Vela X-l, SMC X-l, and Cen X-3 are all V-band data 



Priedhqrskv fc H olt 



1987: 



LMC X-4 is B-band data (jllovaisky et al. 
X-ray eclipse. 



van Paradijs et al. 



19831 ). while the curve for 



19841 ). All data are phased relative to the time of 
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e - 0.174 ± 0.015 e = 




Orbital Phase Orbital Phase 



Fig. 9. — Optical BVI light curves and radial velocity curves for 4U 1538-52. The solid 
lines represent the best fit ELC model, which includes an accretion disk around the neutron 
star. The dotted lines depict the model with the light from the accretion disk subtracted for 
comparison. The left column shows a model with e = 0.174 ± 0.015 and the right column 
shows a model for e = 0. The eccentric orbit model has a lower overall x 2 = 808 (versus 
X 2 = 818 for the circular orbit), but as discussed in Section fl~2l there is no single physical 
solution for an eccentric orbit and a sufficiently long eclipse duration. All data are phased 
relative to the time of X-ray eclipse. 
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Fig. 10. — Final derived neutron star masses for all six systems. The solid circles represent 
masses derived using ELC and the open circles are the analytic solutions for comparison. The 
solid triangle for 4U 1538-52 represents the ELC-derived mass if the system has a circular 
orbit (e = 0), and the open circle for this system also represents a circular orbit (analytic 
solution). The solid circle for this system is the eccentric orbit mass derived using ELC. All 
values are given in Table HJ For consistency, the fit value of j3 in the ELC models was used 
in each analytic solution. We assume /3 = 1 for Her X-l. 
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Table 1. Input Parameters 



System 


P orb (days) 


ax sin i (lt-s) 


e a 


e e (deg) 


K opt (km s x ) 


n 


Ref. 


Vela X-l 


8.964368 ± 0.000040 


113.89 ± 0.13 


0.0898 ± 0.0012 


34.135 ± 0.5 b 


21.7 ± 1.6 


0.67 ±0.04 


1,3,7 


4U 1538-52 


3.728382 ± 0.000011 


56.6±0.7 C 


0.174 ±0.015 c 


28.5 ± 1.5 


20 ±3 


0.91±0.20 d 


2,4,6 


SMC X-l 


3.89229090 ± 0.00000043 


53.4876 ± 0.0004 





28.25 ±2.25 


20.2 ± 1.1 


0.91 ± 0.20 


5 


LMC X-4 


1.40839776 ± 0.00000026 


26.343 ± 0.016 





27 ±2 


35.1 ± 1.5 


0.97 ±0.13 


5 


Cen X-3 


2.08713845 ± 0.00000005 


39.56 ± 0.07 





32.9 ± 1.4 


27.5 ± 2.3 


0.75 ± 0.13 


5 


Her X-l 


1.700167720 ± 0.000000010 


13.1831 ± 0.0003 





24.5 ± 0.5 


90 ± 20 


1.0 ±0.15° 


6 



a Eccentricities less than 0.01 are approximated to be 0. 
b For reasons discussed in Section |4. II we adopt 33 ± 3 from 



van Kerkwiik et al 



(1995bJ) for the analytic case instead. 



c As discussed in Sections 2.1 and 1.2 we consider both an eccentric orbit and the case e = for 4U 1538-52. When e = 0, ax sin i = 54.3±0.6 



i Clark 



2000) 



d We use the same Q here as for SMC X-l since the systems' t> ro t sini values are equal within uncertainty (see Tabled- 
c Here we assume synchronous rotation (f2 = 1). 



References. — (1) Barziv et al. 2001; (2) Clark 2000; (3) Kreykenbohm et al. 2008; (4) Mukherjee et al. 2007; (5) van der Meer et al. 2007 
(and references therein); (6) van Kerkwijk et al. 1995a; (7) Zuiderwijk 1995. 
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Table 2. Additional Input Parameters 



System 


v rot sin i (km s r ) 


uj (degrees) 3. 


Ref. 


Vela X-l 


116 ±6 


332.59 ±0.92 


1,6 


4U 1538-52 


180 ± 20 


198 ± 14 b 


2,3,4 


SMC X-l 


170 ± 30 




5 


LMC X-4 


240 ± 25 




5 


Cen X-3 


200 ± 40 




5 


Her X-l 









a u is the argument of periastron for the companion star 
which is defined only for the systems with nonzero eccen- 
tricity (see Table [1]). 

b This value assumes a constantly changing u over time, 
as discussed in Section 14.21 

References. — (1) Bildsten et al. 1997; (2) Clark 2000; 
(3) Mukherjee et al. 2007; (4) Reynolds et al. 1992; (5) 
van der Meer et al. 2007; (6) van Kerkwijk et al. 1995a. 
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Table 3. Analytic Neutron Star Masses' 



System 


M x (M ) 


% (deg) 


Vela X-l 


1.617 ±0.130 


85.9 ±2.0 


4U 1538-52 b 


0.859 ±0.200 


67.3 ±4.4 


SMC X-l 


1.067 ±0.116 


67.8 ±4.7 


LMC X-4 


1.252 ±0.108 


68.8 ±3.9 


Cen X-3 


1.349 ±0.146 


72.5 ±4.8 


Her X-l 


0.890 ±0.280 


82.0 ±3.7 



a All values listed are from Monte Carlo 
simulations that draw from a random dis- 
tribution of 0.9 < < 1. 

b Values for 4U 1538-52 assume e = 
in order to arrive at a physical solution as 
discussed in Sections 12.11 and 14.21 
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Table 4. Derived System Parameters 



Analytic 3 Numerical 





M x (Mq) 


i (deg) 


M X (Mq) 


i (deg) 


A/opt (Mq) 




r 


P 


Vela X-l 


1.788 ±0.157 


83.6 ±3.1 


1.770 ±0.083 


78.8 ± 1.2 


24.00 ±0.37 


31.82 ± 0.28 


0.99 ±0.01 


1 


4U 1538-52 (ecc) 






0.874 ± 0.073 


68.0 ± 1.4 


20.72 ± 2.27 


15.72 ± 0.52 


0.88 ±0.02 


0.95 


4U 1538-52 (circ) 


1.104 ±0.177 


72.6 ±4.2 


0.996 ±0.101 


76.8 ±6.7 


14.13 ± 2.78 


12.53 ± 2.11 


0.76 ±0.02 


0.88 


SMC X-l 


1.064 ±0.105 


67.8 ±4.2 


1.037 ±0.085 


68.5 ± 5.2 


15.35 ± 1.53 


15.70 ± 1.36 


0.86 ±0.07 


0.95 


LMC X-4 


1.249 ±0.094 


68.8 ±3.3 


1.285 ±0.051 


67.0 ± 1.9 


14.96 ± 0.58 


7.76 ± 0.32 


0.86 ±0.03 


0.95 


Cen X-3 


1.473 ±0.143 


67.5 ±3.2 


1.486 ± 0.082 


66.7 ±2.4 


22.06 ± 1.37 


12.56 ± 0.56 


> 0.96 


1 


Her X-l 


1.036 ±0.311 


80.5 ±3.8 


1.073 ±0.358 


> 85.9 


2.03 ± 0.37 


3.76 ± 0.54 




1 



a For consistency, these analytic values are derived using the /3 values returned by the numerical model rather than a distribution of 
0-9 < /3 < 1 as in Tableland Figured] 

b The sphere-equivalent radius of the companion star (e.g., see Figure 0. 

c The ELC fill factor, defined as the distance from the companion star's center of mass to the point of the star closest to LI. A value 
of / maps directly to a value for /3, which is the Roche lobe filling factor expressed in terms of the sphere-equivalent volume radius. For 
the eccentric systems, both / and f} are defined at periastron. 
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Table 5. Varying K opt and 6 e in SMC X-l 



M x (M ) 


K opt (km s 


e (deg) 


1.091 ±0.062 


18.0 


55.0 


0.989 ±0.056 


18.0 


58.0 


0.970 ±0.092 


18.0 


61.0 


0.904 ±0.050 


18.0 


64.0 


1.103 ±0.061 


20.0 


55.0 


1.099 ±0.070 


20.0 


58.0 


1.002 ±0.050 


20.0 


61.0 


0.944 ±0.055 


20.0 


64.0 


1.203 ±0.064 


22.0 


55.0 


1.117 ± 0.096 


22.0 


58.0 


1.058 ±0.047 


22.0 


61.0 


1.055 ±0.048 


22.0 


64.0 


1.260 ±0.072 


24.0 


55.0 


1.161 ±0.049 


24.0 


58.0 


1.119 ±0.067 


24.0 


61.0 


1.078 ±0.046 


24.0 


64.0 



a These K opt and 6 e values have the same 
uncertainties as the actual K opt and 9 e for 
SMC X-l, which are given in Tabled] 



